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FOREWORD 

This document is the final report describing the results of effort by 
personnel of Lockheed Missiles & Space Company, Inc., Huntsville Research & 
Engineering Center, for the National Aeronautics and Space Administration 
under Contract NAS8-3281, "Fluid Dynamics Numerical Analysis." The con- 
tractual effort described in this document was performed during the year 
from January 1983 to January 1984. The NASA Technical Director for this 
contract is Dr. William Oran, NASA Headquarters, Washington, D.C. 
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INTRODUCTION AND SUMMARY 

The research program described in this final report is a continuation 
of effort initiated in August 1979 to investigate natural convection in 
various materials processing experiment configurations under microgravity 
conditions simulating the orbiting Space Station environment. In the years 
since this research program was initiated, a number of microgravity natural 
convection problems have been investigated, including the upper limit on 
Rayleigh number for linear (low Rayleigh number) theory to be valid, the 
effect of container shape on convection velocities, and the effect of 
changes in the microgravity vector (direction and magnitude) on convection 
velocities. In addition, the development of natural convection was simu- 
lated in the Lal/Kroes experiment. These results are described in annual 
reports for this contract (Refs. 1 through 3). 

The results described in this document are for the period from January 
1983 to January 1984. During this period,, the effort was concerned pri- 
marily with numerically simulating thermocapillary convection in Dr. Robert 
F. Dressier' s air-jet Marangoni flow suppressing technique. Dr. Dressier, 
of George Washington University, is developing this technique under a sepa- 
rate NASA contract. In cooperation with Dr. Dressier, we performed numer- 
ical simulations for flow under various conditions of fluid properties, 
gravity, and temperature gradient. Results presented herein include com- 
puter generated plots of streamlines, velocities, and temperatures through- 
out the contained fluid flow, and plots of Marangoni flow velocities over 
the free surface. 
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DISCUSSION 


Thermocapillary convection was simulated for flow in the idealised two- - 
dimensional container illustrated in the following sketch. The container is 



Idealized Two-Dimensional Square Container 
for Thermocapillary Flow Simulation 

square-shaped with the interior dimensions 0.3 cm x 0.3 cm. The top and 
bottom of the container are perfect thermal conductors held at constant tem- 
peratures, with the top at a temperature, AT, greater than the bottom. One 
side of the container is a perfect thermal insulator, and the other side is 
open, so that the liquid surface on that side is free. For these numerical 
simulations, we assume that the liquid free surface is always flat. We also 
assume that there is no heat transfer across tne liquid free surface. Where 
gravity is present, the direction is always downward. 
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A slc.pllfled analysis was performed to develop closed form algebraic 
equations approximately relating fluid velocities witn the various experi- 
mental parameters. This was done to provide Initial estimates for selecting 
experimental conditions to Le numerically simulated. A brief description 
of the development of these simplified equations is given in the appendix. 

The maximum fluid velocity, v , , at steady stute is given by 

max 


v 

max 


T L 
w 

4p 


( 1 ) 


where •{ is the surface shear caused by surface tension gradient, L is the 

length of the side of the square enclosure and y is viscosity. The surface 

shear t is given by 
w 


a da dT s AT 
T w ~ dT dy L 

where da/dT is the surface tension temperature coefficient and AT is the 
temperature difference across the free surface. 


The equilibrium temperature distribution for the fluid-at-rest consists 
of straight-line isotherms parallel to tne x axis. The distortion, Sy, in 
the isotherms due to fluid motion is given by 


6y _ 

L 


1 

2 

12 + it 


v 

max 


L 


a 


(3) 


where a is the thermal diffusivity. The distortion 6y is the displacement 
of the V = T q isotherm from the y = 0.15 cm line ac the free surface (see 
sketch, next page). 


Some type of silicone oil was originally intended for use in the 
experimental program. Silicone oils are available in a wide range of 
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Schematic Showing Displacement of Isotherm 
Due to Fluid Motion 


viscosity grades. The salient properties of Dow Corning 200 fluid (Ref. 4) 
are listed in Table 1. These listed properties are essentially the same 
regardless of the viscosity grade. The surface tension temperature coeffi- 
cient was provided by Dr. Dressier. 

In order to permit reasonable accuracy in the experimental measure- 
ments, the temperature difference AT across the fluid was required to be at 
least 1 C and the free surface fluid velocity at least 1 mm/sec. The first 
set of conditions that we considered was for AT = 5 C using 1000 centistoke 

viscoscity fluid. The surface shear, t , for this temperature difference 

W 2 

was estimated from Eq. (2) to be about. -1.0 dyne/cm . The maximum fluid 

free surface velocity v was then estimated from Eq. (1) to be about 

max ' * 

-0.07 mm/ sec (flowing downward), considerably less than the 1 mm/sec minimum 
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In our experimental guidelines. The displacement of the Isotherms was esti- 
mated from Eq. (3) to be approximately 8 percent of the distance across the 
free surface. This indicates that the temperature distribution across the 
free surface should not be perturbed greatly by the fluid motion. 

We numerically simulated these conditions using a Navier-Stokes code 
which included surface shear and gravity. The computations were made on the 
CYBER 203 computer at NASA-Langley Research Center. The Prandtl number for 
the 1000 centistoke oil is 8570, which is the ratio of the vorticity 
dlffusivity (kinematic viscosity) to the thermal diffusivity. This high 
ratio of diffusivities indicates that the vorticity (velocity) field will be 
diffused at a much higher rate than the temperature field. Thus, a quasi- 
steady velocity field will be established long before the equilibrium tem- 
perature field is established. Since the temperature field affects the 
velocity field through the surface shear boundary condition, the quasi'* 
steady velocity field slowly changes while the temperature field approaches 
its steady state. This disparity in the two dif f usivities creates a diffi- 
culty in the numerical computations. The integration time step, At, is 
given by 


At 


, . 2 

1 / Ax N 

2 a + v 


(4) 


where Ax is the spacing distance between grid lines, a is thermal dif- 
fusivity and V is vorticity diffusivity (kinematic viscosity). Thus, the 
integration time step is limited by the larger of the two dif f usivities . In 
order to maintain stability, a very small time step is required for the 
coupled vorticity and thermal diffusion problem compared to what would be 
required for the thermal diffusion problem alone. Thus for a fully coupled 
solution, ax^ extremely large number of integration time steps would be 
required to reach the steady state solution. 

This numerical difficulty was dealt witn by resorting to an alternating 
coupled-uncoupled solution technique. A fully coupled solution, using the 
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small time step required for the coupled solution, was used to bring the 
velocity field to a quasi-steady state condition. Then, holding the veloc- 
ity field fixed, the temperature advection-diffusion equation was integrated 
at a much larger time step based on the thermal diffusivity. This procedure 
was alternated until convergence was achieved. 

The results of the numerical integration are shown graphically in Figs. 
1 through 4. These figures present, respectively, streamlines, absolute 
velocity contours, velocity vectors and temperature contours at steady 
state. The maximum computed velocity is about 0.05 mm/sec compared to the 
0,07 mm/sec estimated from Eq, (1). The mid-temperature contour is shown to 
be distorted by about 4 percent of the distance across the free surface due 
to fluid motion, compared to the 8 percent predicted from Eq. (3). 

The good agreement between the numerically computed results and esti- 
mates obtained by simple algebraic equations essentially verifies the accu- 
racy of the computer program. Unfortunately, the computed velocities for 
this first set of conditions (1000 centistolce viscosity grade oil and 5 C 
temperature difference) are well below the 1 mm/ sec experimental guidelines. 

Referring to Eq. (3), we find that, at the minimum desirable fluid 
velocity of 1 mm/sec, the degree of distortion in the temperature contours 
is unacceptably high (greater than 100 percent). The problem lies in the 
extremely low thermal diffusivity of the silicone oil. At this point, we 
have arrived at a dilemma. It appears that we cannot have acceptably nigh 
fluid velocities without producing large distortions in the temperature 
field, using the silicone oil. If the silicone oil is used, then large dis- 
tortions in the temperature field will have to be accepted and dealt with in 
the experimental program. 

Several other sets of candidate experimental conditions were simulated 
to determine what effect the various parameters will have on the flow 
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velocities and temperature field. An extreme example is the case cf the 10 
centistoke viscosity grade fluid and 5 C temperature difference. Computer 
generated plots of streamlines, absolute velocity and temperature contours, 
and velocity vectors are shown in Figs. % through 8. Note the large dis-* 
tortions exhibited in these results compared to the symmetry in the results 
ghown in Figs. 1 through 4. The high fluid velocities in the 10 centistpke 
oil resulted in large distortions in the temperature field. The distortions 
in the temperature field produced increased temperature gradients toward the 
downstream end of the free surface, which produced still higher surface 
shears and, hence, higher velocities. The velocity near the downstream end 
of the free surface, therefore, becomes much higher than the average across 
the surface, apparently creating the vortex shown in the streamline plot in 
Fig. 5. It would obviously be extremely difficult to accurately measure 
velocities in this kind of distribution. 

We next tried 50 centistoke oil with a temperature difference of 2 C. 
These results are shown in Figs. 9, 10, and 11. The streamlines are shown 
in Fig, 9 to be much less distorted than in the 10 centistoke oil. The max- 
imum velocity shown in Fig. 10 is about 1 ram/sec, which corresponds well 
with our experimental guidelines. Here again, however, the higher veloc- 
ities are concentrated near the downstream end of the free surface, and the 
average velocities over the free surface are much lower. This is indicated 
also in the temperature distribution shown in Fig. 11. 

As discussed earlier, the distortions in the velocit / and temperature 
fields are due to the extremely low thermal diffusivity of the silicone oil. 
In order to provide an indication of the behavior of a higher thermal dif- 
fusivity fluid, we increased the thermal diffusivity of the 50 centistoke 
fluid by a factor of 20 and recomputed. The results are shown in Figs. 12 
through 15. Note that the streamlines shown in Fig. 12 and velocity fields 
shown in Figs. 13 and 14 are nearly symmetrical. Also, the amount of dis- 
tortion in the temperature contours is shown to be greatly reduced as com- 
pared to the normal thermal diffusivity silicone oil. 
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The next set of teat conditions that we simulated was for tne 30 ceutl- 
stotce fluid witn a temperature difference of 0.3 C. We anticipated that 
this would produce velocities well below the experimental guidelines, and 
the temperature difference was also somewhat lower than desired from the 
experimental standpoint. The idea behind performing these simulations was 
to observe the amount of velocity and temperature field distortion for 
silicone oil at the extreme low end of the range of conceivably acceptable 
velocities. The results are shown in Figs. 16 through 19, The maximum 
velocity was found to be about 0.1 mm/sec, an order of magnitude below the 
experimental guidelines. The streamlines are shown in Fig. 16 to be 
slightly distorted. The velocity field is shown in Figs. 17 and 18 to 
Increase continuously to a maximum at the downstream end of the free sur- 
face, The temperature field is shown in Fig. 19 to be considerably dis- 
torted even at the small fluid y« k ocities observed for this case. 

The final set of test conditions that we numerically simulated was for 
20 centistoke oil with a temperature difference of 3 C. We anticipated that 
these conditions would result in acceptably high fluid velocities, but would 
display unacceptably high distortions in the temperature and velocity 
fields. These results are shown in Figs. 20 througn 23 for zero gravity and 
in Figs. 24 through 27 with a downward directed gravitational force of 1 g 
imposed. The effect of gravity appears to be negligible. Note the severe 
distortion in the streamline, velocity and temperature fields. Surface 
velocities are essentially the same for both gravity conditions and are 
plotted in Fig. 28 as one curve. The surface velocities are shown to be 
relatively uniform at about 1 mm/sec over a large portion of the free sur- 
face. An extreme value of about 5 mra/sec is reached in a sharp peak near 
the downstream end. 

Again, 'as an indication of what might be expected from a higher thermal 
diffusivlty fluid, we increased the thermal dlffusivity by a factor of 100, 
and repeated the computations for the 20 centistoke fluid with 5 C tempera- 
ture difference. Computer generated plots are shown in Figs. 29 through 32. 
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These plots all show very little distortions in the velocity and temperature 
fields. The plot of surface velocity in Fig, 33 shows a smooth parabolic" 
like curve with o muxlmum velocity of about 2 mm/sec in the center of the 
free surface. Increasing the thermal diffusivity by a factor of 100 places 
it in the range of liquid metalB, e.g. , melted tin. 
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CONCLUSIONS 

The thermocapillary convection experimental configuration was success- 
fully numerically simulated over a wide range of possible test conditions. 
Some difficulty arose as a result of the extremely low thermal diffusivity 
of the Bilicone oil. Fx.at, the wide disparity between the vorticlty and 
thermal diffusion rates resulted in extremely long convergence times in the 
numerical simulation. This difficulty was overcome by the use of alterna- 
ting coupled and uncoupled numerical integrations of the mass, momentum and 
energy conservation equations. Second, the low thermal diffusivity was 
shown to result in undesirable distortions in the velocity and temperature 
fields. If these undesirable distortions are determined to be unacceptable, 
a higher thermal diffusivity fluid will be required. 
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Table 1 PROPERTIES OF SILICONE OIL 


Density, p 

Thermal Conductivity, k 

Specific Heat, C 

Volumetric Coefficient of 
Thermal Expansion, 8 t 

Surface Tension, a 

Surface Tension Temperature 
Coefficient, da/dT 

Kinematic Viscosity, V 


0.971 gm/cm^ 

0.00038 cal/cm-sec-C 
0,335 cal/gm-C 
0.00096/C 

21. 2 dyne /cm 
-0.058 dyne/cm-C 

Available in Wide Range 
of Viscosity Grades 
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Fig. 33 Velocity Along Free Surface for 20 Centistoke 
Silicone Oil with Temperature Difference of 
5 C and Thermal Diffusivity Increased by 
Factor of 100 
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FOR ESTIMATING MARANGONI 
CONVECTION IN SQUARE ENCLOSURES 
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Appendix 


We assume simple algebraic forms for the velocity components: 

u = U(x) sin(2TT — — j~~) 
v = V(x) cos(lT - — 

L 


(A.l) 
(A. 2) 


These equations provide the general circulatory flow pattern that we expect 
from intuition. The conservation of momentum in the y direction is ex- 
pressed by the following differential equation: 


(puv - T) + ~ (pv 2 + p) » 0 


(A. 3) 


After dropping insignificant terms, substituting x = y 3v/3x, and using the 
form of v given by Eq. (A. 2), this becomes: 


y 



cos(tt 



a p 
sy 


(A. 4) 


Integrating this equation with respect to x along y ■ L/2, with the 
boundary conditions, v = 0 at x = 0 and x w = y 3v/3x at x = L, yields: 

v ■ ff {T w - §7 < L - y» * < A - 5 > 


A-l 
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This equation automatically satisfies force balances on the overall system 
Ln Lho y-diruct ion, but it contains the unknown 0 p/0 y . This is removed by 
satisfying a mass balance along y * L/2: 


This results in 



9p 3 

* 'T 

3y 2L w 


(A. 6) 


(A.7) 


Combining Eqs. (A. 2), (A. 5), and (A.7) results in the expression for the 
v-component of velocity: 


v 


T L 


w 


2d 


(| 5- 1) £cos(„ 



(A. 8) 


with the maximum v-component velocity, v max > being at the free surface (x 
= L) and given by 


v 

max 


T w L 

Ay 


(A. 9) 


Note that, for x w > 0, v is negative for x < 2/3 L, positive for x > 

2/3 L and, of course, zero at x = 2/3 L. The general flow, therefore, 
obviously is counterclockwise around the point x « 2/3 L and y - L/2. Mass 
conservation requires that: 


/ 

2/3 L 


v dx 


y=L/2 


/ 

1/2 L 


u dy 


x=2/3 L 


(A. 10) 


A- 2 


LOCKHEED-HUNTSVILLE RESEARCH & ENGINEERING CENTER 


ORIGINAL PAGE fS 
OP POOR QUALITY 


LMSC-HRKC TR D951333 


We assume U(x) in Eq. (A.l) to be of the form 


A(1 - £> (i) 


(A. 11) 


This provides an extreme value for v at x ** 2/3 L and zeros at x ° 0 and L. 
Carrying out the integrations in Eq. (A. 10) yields 


A 


T L 
2L 

4 U 


(A. 12) 


Combining Eqs. (A.l), (A. 11), and (A. 12) yields the expression for the 
u-component of velocity: 


u 



T L 
w 


(1 




sin (27T 


with the maximum u-component velocity, u , being at x » 2/3 L and 
LI max 

y = j (1 + j) and given by 

7T T W L _ A7T_ 

u max 27 y 27 v max 


The temperature field may be represented by 


(A. 13) 


(A. 14) 


T(x,y) = T(y) + T’(x,y) (A. 15) 

where T is the temperature field for the fluid at rest and T f is the per- 
turbation due to fluid motion. 

The fluid-at-rest temperature distribution is given by 


t = T 0 + (y " L/2 > 


(A. 16) 
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OF POOR 

where T Q Is the mid-point or average fluid temperature and AT is the tem- 
perature difference across the fluid. 

We assume T* to be of the form 

2 

T’ * B[<1- !£)<£) - 1] cos ( 7T (A. 17) 

This provides for the proper boundary conditions BT*/3x * 0 at x « 0 and L, 
and T' « 0 at y - 0 and L. It also results in 

L L 

/ T' dx dy » 0 (A. 18) 

o o 



This maintains the original average fluid temperature T . We find the 
factor B by requiring the satisfaction of the advection-dif fusion equation 
at the free surface mid-point (x = L and y ■ L/2): 


u|I + V 

Bx 


3T 

By 


2 2 

_ „ JT , 9^T s 

— 01 ( „ + o ) 


Bx* 


By 


(A. 19) 


After inserting the values of u, v and the T derivatives at x = L and y = 
L/2, we find 


B 


v L 
max 


2(1 + 7T / 12 ) 


AT 


The expression for the temperature field then becomes 
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(A. 21) 


By setting the temperature T in Eq. (A. 21) equal to the fluld-at-rest mid 
point temperature T , the displacement, 6y » y - W2, of the mid-isotherm 
at the free surface (x » L) is found to be: 
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